library(MatchIt)
library(foreign)
library(car)
library(descr)
library(ggplot2)
library(plotrix)
library(plyr)
library(dplyr)
library(haven)
library(graphics)
library(effects)
library(data.table)

beat_data <- read_dta("beat_counts.dta") %>% na.omit()

# Supplemental analysis (Table A13-A16 and Figure A2 of the Appendix)

# convict is the treatment; 14.44 is one standard deviation above the mean
beat_data<- beat_data %>% mutate(con_treat=ifelse(convict_pk>14.44,1,0))

m.out1 <- matchit(con_treat ~ pct_pov+pct_unemp+pct_owner+pct_black+pct_latino+pct_young+
                    pct_old+pct_college,
                  data = beat_data, method = "nearest", distance = "logit")
summary(m.out1)
m.data1 <- match.data(m.out1,distance ="pscore")
# difference in means test
t.test(m.data1$to12[m.data1$con_treat==1],m.data1$to12[m.data1$con_treat==0],paired=TRUE)
t.test(m.data1$to14[m.data1$con_treat==1],m.data1$to14[m.data1$con_treat==0],paired=TRUE)
t.test(m.data1$to15[m.data1$con_treat==1],m.data1$to15[m.data1$con_treat==0],paired=TRUE)
t.test(m.data1$mtg_p1k[m.data1$con_treat==1],m.data1$mtg_p1k[m.data1$con_treat==0],paired=TRUE)
t.test(m.data1$call311_p1k[m.data1$con_treat==1],m.data1$call311_p1k[m.data1$con_treat==0],paired=TRUE)
# interaction between treatment and presence of CBOs 
# (ie: how does CBO density matter for participation outcomes among treated and non-treated beats?)
m1 <- lm(to12 ~ con_treat*lcbo, data = m.data1)
m2 <- lm(to14 ~ con_treat*lcbo , data = m.data1)
m3 <- lm(to15 ~ con_treat*lcbo , data = m.data1)
m4 <- lm(mtg_p1k ~ con_treat*lcbo, data = m.data1)
m5 <- lm(call311_p1k ~ con_treat*lcbo, data = m.data1)


stargazer(m2,m3,m4,m5, no.space=TRUE,type="text",
          single.row = TRUE,
          covariate.labels=c("Log(Convictions per 1000 Pop)",
                             "log(CSO per 1000 Pop)",
                             "Convict*CSOs",
                             "Constant"),
          dep.var.labels.include = FALSE,
          model.numbers=FALSE,
          star.char=c("*","**","***"),
          star.cutoffs = c(0.05,0.01,0.001),
          column.labels = c("2014 Turnout","2015 Turnout",
                            "Meeting Attendence per 1000 Pop","311 Calls Per 1000 Pop"),
          dep.var.labels   = "")

# Visualize
#TO 14
coefs = as.data.frame(summary(m2)$coefficients[-1,1:2])
names(coefs)[2]="se"
coefs$vars = rownames(coefs)

variable<-c("01","02","03")
d<-cbind(variable,coefs)
p1<-ggplot(data = d, aes(x = variable, y = Estimate, 
                         ymin=(Estimate - 1.96*se), ymax = (Estimate + 1.96*se))) +
  geom_pointrange(size=.5)+
  theme_bw(base_size=16,base_family="Times")+
  scale_x_discrete(breaks=c("01","02","03"),labels=c("High Conviction","Log(CSO)","Convict*CSO"))+
  geom_hline(yintercept = 0, lty = 2)+
  theme(plot.margin=unit(c(1,2,1,1),"cm"))+
  labs(y="Marginal Effect\n",x="\n Model: 2014 Turnout",title=c("", cex=1))

ggsave(p1, filename = paste("coeff_14TO_matched", ".pdf", sep = ""),width = 7, height = 5)

# TO 15
coefs = as.data.frame(summary(m3)$coefficients[-1,1:2])
names(coefs)[2]="se"
coefs$vars = rownames(coefs)

variable<-c("01","02","03")
d<-cbind(variable,coefs)
p1<-ggplot(data = d, aes(x = variable, y = Estimate, 
                         ymin=(Estimate - 1.96*se), ymax = (Estimate + 1.96*se))) +
  geom_pointrange(size=.5)+
  theme_bw(base_size=16,base_family="Times")+
  scale_x_discrete(breaks=c("01","02","03"),labels=c("High Conviction","Log(CSO)","Convict*CSO"))+
  geom_hline(yintercept = 0, lty = 2)+
  theme(plot.margin=unit(c(1,2,1,1),"cm"))+
  labs(y="Marginal Effect\n",x="\n Model: 2015 Turnout",title=c("", cex=1))


ggsave(p1, filename = paste("coeff_15TO_matched", ".pdf", sep = ""),width = 7, height = 5)


#meetings
coefs = as.data.frame(summary(m4)$coefficients[-1,1:2])
names(coefs)[2]="se"
coefs$vars = rownames(coefs)

variable<-c("01","02","03")
d<-cbind(variable,coefs)
p1<-ggplot(data = d, aes(x = variable, y = Estimate, 
                         ymin=(Estimate - 1.96*se), ymax = (Estimate + 1.96*se))) +
  geom_pointrange(size=.5)+
  theme_bw(base_size=16,base_family="Times")+
  scale_x_discrete(breaks=c("01","02","03"),labels=c("High Conviction","Log(CSO)","Convict*CSO"))+
  geom_hline(yintercept = 0, lty = 2)+
  theme(plot.margin=unit(c(1,2,1,1),"cm"))+
  labs(y="Marginal Effect \n",x="\n Model: Attendance at Police Beat Meetings",title=c("", cex=1))

ggsave(p1, filename = paste("coeff_meetings_matched", ".pdf", sep = ""),width = 7, height = 5)

# 311 calls
coefs = as.data.frame(summary(m5)$coefficients[-1,1:2])
names(coefs)[2]="se"
coefs$vars = rownames(coefs)

variable<-c("01","02","03")
d<-cbind(variable,coefs)
p1<-ggplot(data = d, aes(x = variable, y = Estimate, 
                         ymin=(Estimate - 1.96*se), ymax = (Estimate + 1.96*se))) +
  geom_pointrange(size=.5)+
  theme_bw(base_size=16,base_family="Times")+
  scale_x_discrete(breaks=c("01","02","03"),labels=c("High Conviction","Log(CSO)","Convict*CSO"))+
  geom_hline(yintercept = 0, lty = 2)+
  theme(plot.margin=unit(c(1,2,1,1),"cm"))+
  labs(y="Marginal Effect\n",x="\n Model: 311 calls per 1k pop",title=c("", cex=1))

ggsave(p1, filename = paste("coeff_311calls_matched", ".pdf", sep = ""),width = 7, height = 5)

